source("include.R")
❤ Loaded agrimonia dataset. Available as df_agri.
❤ Converted df_agri$Time to Date variable type, year-month-day.
❤ Created DATE_FORMAT variable for date comparisons.
❤ Loaded stations split dataset. Available as df_stat.

Maps plots

library(maps)
Avvertimento: il pacchetto ‘maps’ è stato creato con R versione 4.2.3
library(ggplot2)

Caricamento pacchetto: ‘ggplot2’

Il seguente oggetto è mascherato da ‘package:crayon’:

    %+%
library(sf)
Avvertimento: il pacchetto ‘sf’ è stato creato con R versione 4.2.3Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE

Idea: plot the values for pm10 each day, like with “bubbles” of ggplot, so to get 2192 plots. And then combine all them in a single gif, which shows the evolution over time of those pm10 values.

https://askubuntu.com/questions/648244/how-do-i-create-an-animated-gif-from-still-images-preferably-with-the-command-l

So idea: plot bubbles with size according to the pm10 value. And maybe if we want use also a different color (like to separate high/medium/low values).

regioni_italiane <- st_read("italia/gadm40_ITA_1.shp") 
Reading layer `gadm40_ITA_1' from data source 
  `C:\Users\feder\Desktop\Uni magistrale\Bayesian statistics\progetto-bayesian\src\italia\gadm40_ITA_1.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 20 features and 11 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 6.630879 ymin: 35.49292 xmax: 18.52069 ymax: 47.09096
Geodetic CRS:  WGS 84
sites <- data.frame(
    longitude = unique(df_agri$Longitude), 
    latitude = unique(df_agri$Latitude))

pad = 0.2*c(-1,1)
mappa_completa <- ggplot() +
  geom_sf(data = regioni_italiane, fill = "lightblue", color = "black", size = 0.5)+
      coord_sf(xlim = range(sites$longitude)+pad, 
             ylim = range(sites$latitude)+pad, expand = FALSE)+
  # geom_sf(data = world) +
  geom_point(data = sites, aes(x = longitude, y = latitude), size = 1, shape = 23, fill = "darkred")
 
print(mappa_completa)

sites <- data.frame(
    longitude = unique(df_agri$Longitude), 
    latitude = unique(df_agri$Latitude))

for (data in seq.Date(from=as.Date("2016-01-01",DATE_FORMAT),
               to=as.Date("2016-01-11",DATE_FORMAT),by="days")){
        
    dt = as.Date(data,origin="1970-01-01")
    df_date = df_agri[which(df_agri$Time == dt),]
    print(dt)
    
    pad = 0.2*c(-1,1)
    mappa_completa <- ggplot() +
      geom_sf(data = regioni_italiane, fill = "lightblue", color = "black", size = 0.5)+
          coord_sf(xlim = range(sites$longitude)+pad, 
                 ylim = range(sites$latitude)+pad, expand = FALSE)+
      # geom_sf(data = world) +
      geom_point(data = sites, aes(x = longitude, y = latitude), size = 1, shape = 23, fill = "darkred")+
        geom_point(data=df_date,aes(x=Longitude,y=Latitude),size=df_date$AQ_pm10/10,alpha=0.5)+
        ggtitle(dt)
    
    print(mappa_completa)
    ggsave(paste0("./maps_images/",dt,".jpeg"),plot=mappa_completa)
}
[1] "2016-01-01"
[1] "2016-01-02"
[1] "2016-01-03"
[1] "2016-01-04"
[1] "2016-01-05"
[1] "2016-01-06"
[1] "2016-01-07"
[1] "2016-01-08"
[1] "2016-01-09"
[1] "2016-01-10"
[1] "2016-01-11"

Trend plots

stations = unique(df_agri$IDStations)
cols = colora(6,97,show=1)


for (st in stations[1:4]){
    df_st = df_stat[[st]]
    data2016 = subset(df_st,Time>="2016-01-01" & Time<="2016-12-31")
    data2017 = subset(df_st,Time>="2017-01-01" & Time<="2017-12-31")
    data2018 = subset(df_st,Time>="2018-01-01" & Time<="2018-12-31")
    data2019 = subset(df_st,Time>="2019-01-01" & Time<="2019-12-31")
    data2020 = subset(df_st,Time>="2020-01-01" & Time<="2020-12-31")
    data2021 = subset(df_st,Time>="2021-01-01" & Time<="2021-12-31")
    
    
    plot(as.Date(format(data2016$Time,"%d %m"),"%d %m"),data2016$AQ_pm10 ,type="l",col=cols[1],
         xlab="time",ylab="PM_10 values",main=paste0("Station ",st,", all years"),
         ylim=range(na.omit(df_st$AQ_pm10)))
    lines(as.Date(format(data2017$Time,"%d %m"),"%d %m"),data2017$AQ_pm10,type="l",col=cols[2])
    lines(as.Date(format(data2018$Time,"%d %m"),"%d %m"),data2018$AQ_pm10,type="l",col=cols[3])
    lines(as.Date(format(data2019$Time,"%d %m"),"%d %m"),data2019$AQ_pm10,type="l",col=cols[4])
    lines(as.Date(format(data2020$Time,"%d %m"),"%d %m"),data2020$AQ_pm10,type="l",col=cols[5])
    lines(as.Date(format(data2021$Time,"%d %m"),"%d %m"),data2021$AQ_pm10,type="l",col=cols[6])
    legend('topright', levels(as.factor(2016:2021)), fill=cols, bty='n', cex = 0.6)

}

cols = colora(10,show=0)
for (year in 2016:2021){
    for (i in 1:length(stations)){
        st = stations[i]
        df_st = df_stat[[st]]
        df_st = df_st[df_st$Time>=as.Date(paste0(year,"-01-01"),"%Y-%m-%d") &
                      df_st$Time<=as.Date(paste0(year,"-12-31"),"%Y-%m-%d"),]
        if(st=="1264"){
            plot(df_st$Time,df_st$AQ_pm10,type="l",col=cols[i%%10],
                 ylim=range(na.omit(df_agri$AQ_pm10)),main=paste0(year," all the stations"),
                 xlab="months",ylab="PM_10 values")
        } else{
            lines(df_st$Time,df_st$AQ_pm10,type="l",col=cols[i%%10])
        }
    }
}

LS0tDQp0aXRsZTogIlIgTm90ZWJvb2siDQpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sNCi0tLQ0KDQpgYGB7cn0NCnNvdXJjZSgiaW5jbHVkZS5SIikNCmBgYA0KDQojIE1hcHMgcGxvdHMNCmBgYHtyfQ0KbGlicmFyeShtYXBzKQ0KbGlicmFyeShnZ3Bsb3QyKQ0KbGlicmFyeShzZikNCmBgYA0KDQpJZGVhOiBwbG90IHRoZSB2YWx1ZXMgZm9yIHBtMTAgZWFjaCBkYXksIGxpa2Ugd2l0aCAiYnViYmxlcyIgb2YgZ2dwbG90LCBzbyB0byBnZXQgMjE5MiBwbG90cy4gQW5kIHRoZW4gY29tYmluZSBhbGwgdGhlbSBpbiBhIHNpbmdsZSBnaWYsIHdoaWNoIHNob3dzIHRoZSBldm9sdXRpb24gb3ZlciB0aW1lIG9mIHRob3NlIHBtMTAgdmFsdWVzLg0KDQpodHRwczovL2Fza3VidW50dS5jb20vcXVlc3Rpb25zLzY0ODI0NC9ob3ctZG8taS1jcmVhdGUtYW4tYW5pbWF0ZWQtZ2lmLWZyb20tc3RpbGwtaW1hZ2VzLXByZWZlcmFibHktd2l0aC10aGUtY29tbWFuZC1sDQoNClNvIGlkZWE6IHBsb3QgYnViYmxlcyB3aXRoIHNpemUgYWNjb3JkaW5nIHRvIHRoZSBwbTEwIHZhbHVlLg0KQW5kIG1heWJlIGlmIHdlIHdhbnQgdXNlIGFsc28gYSBkaWZmZXJlbnQgY29sb3IgKGxpa2UgdG8gc2VwYXJhdGUgaGlnaC9tZWRpdW0vbG93IHZhbHVlcykuDQoNCmBgYHtyfQ0KcmVnaW9uaV9pdGFsaWFuZSA8LSBzdF9yZWFkKCJpdGFsaWEvZ2FkbTQwX0lUQV8xLnNocCIpIA0KDQpzaXRlcyA8LSBkYXRhLmZyYW1lKA0KCWxvbmdpdHVkZSA9IHVuaXF1ZShkZl9hZ3JpJExvbmdpdHVkZSksIA0KCWxhdGl0dWRlID0gdW5pcXVlKGRmX2FncmkkTGF0aXR1ZGUpKQ0KDQpwYWQgPSAwLjIqYygtMSwxKQ0KbWFwcGFfY29tcGxldGEgPC0gZ2dwbG90KCkgKw0KICBnZW9tX3NmKGRhdGEgPSByZWdpb25pX2l0YWxpYW5lLCBmaWxsID0gImxpZ2h0Ymx1ZSIsIGNvbG9yID0gImJsYWNrIiwgc2l6ZSA9IDAuNSkrDQoJICBjb29yZF9zZih4bGltID0gcmFuZ2Uoc2l0ZXMkbG9uZ2l0dWRlKStwYWQsIA0KCSAgCQkgeWxpbSA9IHJhbmdlKHNpdGVzJGxhdGl0dWRlKStwYWQsIGV4cGFuZCA9IEZBTFNFKSsNCiAgIyBnZW9tX3NmKGRhdGEgPSB3b3JsZCkgKw0KICBnZW9tX3BvaW50KGRhdGEgPSBzaXRlcywgYWVzKHggPSBsb25naXR1ZGUsIHkgPSBsYXRpdHVkZSksIHNpemUgPSAxLCBzaGFwZSA9IDIzLCBmaWxsID0gImRhcmtyZWQiKQ0KIA0KcHJpbnQobWFwcGFfY29tcGxldGEpDQpgYGANCg0KDQpgYGB7cn0NCnNpdGVzIDwtIGRhdGEuZnJhbWUoDQoJbG9uZ2l0dWRlID0gdW5pcXVlKGRmX2FncmkkTG9uZ2l0dWRlKSwgDQoJbGF0aXR1ZGUgPSB1bmlxdWUoZGZfYWdyaSRMYXRpdHVkZSkpDQoNCmZvciAoZGF0YSBpbiBzZXEuRGF0ZShmcm9tPWFzLkRhdGUoIjIwMTYtMDEtMDEiLERBVEVfRk9STUFUKSwNCgkJCSAgIHRvPWFzLkRhdGUoIjIwMTYtMDEtMTEiLERBVEVfRk9STUFUKSxieT0iZGF5cyIpKXsNCgkJDQoJZHQgPSBhcy5EYXRlKGRhdGEsb3JpZ2luPSIxOTcwLTAxLTAxIikNCglkZl9kYXRlID0gZGZfYWdyaVt3aGljaChkZl9hZ3JpJFRpbWUgPT0gZHQpLF0NCglwcmludChkdCkNCgkNCglwYWQgPSAwLjIqYygtMSwxKQ0KCW1hcHBhX2NvbXBsZXRhIDwtIGdncGxvdCgpICsNCgkgIGdlb21fc2YoZGF0YSA9IHJlZ2lvbmlfaXRhbGlhbmUsIGZpbGwgPSAibGlnaHRibHVlIiwgY29sb3IgPSAiYmxhY2siLCBzaXplID0gMC41KSsNCgkJICBjb29yZF9zZih4bGltID0gcmFuZ2Uoc2l0ZXMkbG9uZ2l0dWRlKStwYWQsIA0KCQkgIAkJIHlsaW0gPSByYW5nZShzaXRlcyRsYXRpdHVkZSkrcGFkLCBleHBhbmQgPSBGQUxTRSkrDQoJICAjIGdlb21fc2YoZGF0YSA9IHdvcmxkKSArDQoJICBnZW9tX3BvaW50KGRhdGEgPSBzaXRlcywgYWVzKHggPSBsb25naXR1ZGUsIHkgPSBsYXRpdHVkZSksIHNpemUgPSAxLCBzaGFwZSA9IDIzLCBmaWxsID0gImRhcmtyZWQiKSsNCgkJZ2VvbV9wb2ludChkYXRhPWRmX2RhdGUsYWVzKHg9TG9uZ2l0dWRlLHk9TGF0aXR1ZGUpLHNpemU9ZGZfZGF0ZSRBUV9wbTEwLzEwLGFscGhhPTAuNSkrDQoJCWdndGl0bGUoZHQpDQoJDQoJcHJpbnQobWFwcGFfY29tcGxldGEpDQoJZ2dzYXZlKHBhc3RlMCgiLi9tYXBzX2ltYWdlcy8iLGR0LCIuanBlZyIpLHBsb3Q9bWFwcGFfY29tcGxldGEpDQp9DQpgYGANCg0KDQoNCiMgVHJlbmQgcGxvdHMNCmBgYHtyfQ0Kc3RhdGlvbnMgPSB1bmlxdWUoZGZfYWdyaSRJRFN0YXRpb25zKQ0KY29scyA9IGNvbG9yYSg2LDk3LHNob3c9MSkNCg0KZm9yIChzdCBpbiBzdGF0aW9uc1sxOjRdKXsNCglkZl9zdCA9IGRmX3N0YXRbW3N0XV0NCglkYXRhMjAxNiA9IHN1YnNldChkZl9zdCxUaW1lPj0iMjAxNi0wMS0wMSIgJiBUaW1lPD0iMjAxNi0xMi0zMSIpDQoJZGF0YTIwMTcgPSBzdWJzZXQoZGZfc3QsVGltZT49IjIwMTctMDEtMDEiICYgVGltZTw9IjIwMTctMTItMzEiKQ0KCWRhdGEyMDE4ID0gc3Vic2V0KGRmX3N0LFRpbWU+PSIyMDE4LTAxLTAxIiAmIFRpbWU8PSIyMDE4LTEyLTMxIikNCglkYXRhMjAxOSA9IHN1YnNldChkZl9zdCxUaW1lPj0iMjAxOS0wMS0wMSIgJiBUaW1lPD0iMjAxOS0xMi0zMSIpDQoJZGF0YTIwMjAgPSBzdWJzZXQoZGZfc3QsVGltZT49IjIwMjAtMDEtMDEiICYgVGltZTw9IjIwMjAtMTItMzEiKQ0KCWRhdGEyMDIxID0gc3Vic2V0KGRmX3N0LFRpbWU+PSIyMDIxLTAxLTAxIiAmIFRpbWU8PSIyMDIxLTEyLTMxIikNCgkNCgkNCglwbG90KGFzLkRhdGUoZm9ybWF0KGRhdGEyMDE2JFRpbWUsIiVkICVtIiksIiVkICVtIiksZGF0YTIwMTYkQVFfcG0xMCAsdHlwZT0ibCIsY29sPWNvbHNbMV0sDQoJCSB4bGFiPSJ0aW1lIix5bGFiPSJQTV8xMCB2YWx1ZXMiLG1haW49cGFzdGUwKCJTdGF0aW9uICIsc3QsIiwgYWxsIHllYXJzIiksDQoJCSB5bGltPXJhbmdlKG5hLm9taXQoZGZfc3QkQVFfcG0xMCkpKQ0KCWxpbmVzKGFzLkRhdGUoZm9ybWF0KGRhdGEyMDE3JFRpbWUsIiVkICVtIiksIiVkICVtIiksZGF0YTIwMTckQVFfcG0xMCx0eXBlPSJsIixjb2w9Y29sc1syXSkNCglsaW5lcyhhcy5EYXRlKGZvcm1hdChkYXRhMjAxOCRUaW1lLCIlZCAlbSIpLCIlZCAlbSIpLGRhdGEyMDE4JEFRX3BtMTAsdHlwZT0ibCIsY29sPWNvbHNbM10pDQoJbGluZXMoYXMuRGF0ZShmb3JtYXQoZGF0YTIwMTkkVGltZSwiJWQgJW0iKSwiJWQgJW0iKSxkYXRhMjAxOSRBUV9wbTEwLHR5cGU9ImwiLGNvbD1jb2xzWzRdKQ0KCWxpbmVzKGFzLkRhdGUoZm9ybWF0KGRhdGEyMDIwJFRpbWUsIiVkICVtIiksIiVkICVtIiksZGF0YTIwMjAkQVFfcG0xMCx0eXBlPSJsIixjb2w9Y29sc1s1XSkNCglsaW5lcyhhcy5EYXRlKGZvcm1hdChkYXRhMjAyMSRUaW1lLCIlZCAlbSIpLCIlZCAlbSIpLGRhdGEyMDIxJEFRX3BtMTAsdHlwZT0ibCIsY29sPWNvbHNbNl0pDQoJbGVnZW5kKCd0b3ByaWdodCcsIGxldmVscyhhcy5mYWN0b3IoMjAxNjoyMDIxKSksIGZpbGw9Y29scywgYnR5PSduJywgY2V4ID0gMC42KQ0KDQp9DQpgYGANCg0KDQpgYGB7cn0NCmNvbHMgPSBjb2xvcmEoMTAsc2hvdz0wKQ0KZm9yICh5ZWFyIGluIDIwMTY6MjAyMSl7DQoJZm9yIChpIGluIDE6bGVuZ3RoKHN0YXRpb25zKSl7DQoJCXN0ID0gc3RhdGlvbnNbaV0NCgkJZGZfc3QgPSBkZl9zdGF0W1tzdF1dDQoJCWRmX3N0ID0gZGZfc3RbZGZfc3QkVGltZT49YXMuRGF0ZShwYXN0ZTAoeWVhciwiLTAxLTAxIiksIiVZLSVtLSVkIikgJg0KCQkJCQkgIGRmX3N0JFRpbWU8PWFzLkRhdGUocGFzdGUwKHllYXIsIi0xMi0zMSIpLCIlWS0lbS0lZCIpLF0NCgkJaWYoc3Q9PSIxMjY0Iil7DQoJCQlwbG90KGRmX3N0JFRpbWUsZGZfc3QkQVFfcG0xMCx0eXBlPSJsIixjb2w9Y29sc1tpJSUxMF0sDQoJCQkJIHlsaW09cmFuZ2UobmEub21pdChkZl9hZ3JpJEFRX3BtMTApKSxtYWluPXBhc3RlMCh5ZWFyLCIgYWxsIHRoZSBzdGF0aW9ucyIpLA0KCQkJCSB4bGFiPSJtb250aHMiLHlsYWI9IlBNXzEwIHZhbHVlcyIpDQoJCX0gZWxzZXsNCgkJCWxpbmVzKGRmX3N0JFRpbWUsZGZfc3QkQVFfcG0xMCx0eXBlPSJsIixjb2w9Y29sc1tpJSUxMF0pDQoJCX0NCgl9DQp9DQpgYGANCg0KDQoNCg0KDQoNCg0KDQoNCg0KDQoNCg==